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!  IiVIRUlnif  !  RiM 


The  general  p'jiri*  U'lii  mat  describe  the  Mow  c.t  the  ocesn  are  well 
known,  but  for  m a  n  .  c  a  s  e  -  of  interest,  analytical  solutions  are  not 
possible  because  of  irregular  boundaries  or  significant  non-linearities 
in  the  equations.  It  is  often  possible  to  obtain  useful  approximate 
solutions  using  numerical  techniques.  However,  for  reasonably  detailed 
numerical  Models,  very  large  amounts  of  computes  tine  are  required, 
and  M  becomes  imperative  to  seek  out  the  Most  efficient  numerical 
algorithms.  In  the  case  where  forecasts  are  being  made  on  a  real  tine 
basis,  the  use  of  non-optimum  schemes  may  result  in  unacceptable 
dr  lavs . 


(bis  report  describes  a  class  of  techniques  for  the  efficient  treatment 
of  bodies  of  water  with  irregular  boundaries.  Realistic  Models  of  the 
ocean  state  must  include  adequate  treatment  of  the  boundaries,  the 
irregular  bottom  topography,  the  shape  of  the  coastline  and  the 
complicated  geometry  of  groups  of  islands. 

Traditional  numerical  methods  using  a  rectangular  finite  difference  grid 

are  suitable  for  exploratory  surveys  of  idealised  problems  where  the  aim 

is  to  see  why  a  particular  system  responds  as  it  does,  and  wl.at  effects 

variations  in  the  parameters  of  the  problem  and  the  forcing  terms  have 

on  tne  outcome. jr  In  order  to  use  these  models  for  complicated 
1  \ 

Metrics.  it  is  \  necessary  to  use  a  very  fine-  grid,  and  a 


correspondi  ng  increase  in  conputer  time,  or  suffer  from  reduced 

•u:  r  i:r  ac  v  . 

'he  methods  described  here  emplov  an  irregular  triangular  grid,  instead 
of  a  regular  rectangular  grid,  so  it  is  possible  to  fit  empirical 
boundaries  with  high  precision  without  using  prohibitively  fine  gnus 
throughout  the  domain  of  the  solution.  The  penalty  for  this  is  a  slight 
increase  in  the  complexity  of  the  algoi ithm. 

Triangular  grids  have  been  extensively  used  for  elastic  and  plastic  flow 
problems  with  the  finite  element  method  (see  for  example,  Strang  and 
Fix,  19, ’3).  These  grids  have  also  been  used  in  Lagrangian  formulations 
of  the  equations  of  hydrodynamics  by  Crowley  (1971).  Boris  et  al  (1975) 
and  Fritts  (1976),  but  the  Lagrangian  method  is  most  suitable  for  flows 
where  the  total  deformation  of  the  fluid  is  small.  Eulenan 
calculations  have  been  performed  by  Sadourny  et  al  (1968),  Uilliamson 
(1968)  and  Thacker  (1977). 


This  report  describes  t-o  tti  a  Methodology  for  constructing  and  using 
tr.eqular  tnanqular  qnds  for  solving  ocean  flow  problens,  and  a  set 
o’'  computer  progransf or  mplenenting,  testing  and  evaluating  the 
nethods.  The  tonplete  set  of  subroutines  is  given  later  as  an  appendix. 
Certain  details  of  the  basic  algorithm  are  best  understood  bv  reference 
to  the  code  I’selt.  The  code  has  been  designed  to  be  highly  nodular, 
so  that  sone  effort  «ust  be  expended  in  describing  the  interface  to  each 
nodule. 

To  n.i.inise  the  flexibility  of  the  code,  it  has  been  written  in  a 
superset  of  FORTRAN  which  includes  nacre  expansions.  By  the  use  of  the 
fiie  inclusion  nacro,  SINLLUDE,  connon  block  Maintenance  is  greatly 
facilitated,  while  the  bulk  of  the  source  codp  is  reduced.  Many 
conpilers  support  sone  forn  of  file  inclusion.  The  other  use  of  Macros 
in  the  code  is  to  support  paraneters.  Paraneters  are  used  for  values 
which  mist  take  the  torn  of  constants,  for  exanple  to  dimension  an 
array.  The  values  nay  need  to  be  changed  to  produce  different  versions 
of  the  rode,  for  exanple  to  be  able  to  change  the  nunber  of  grid  points 
or  to  handle  different  hardware  configurations.  The  easiest  way  to  be 
able  to  enploy  Macros  is  to  use  a  sinple  pre-processor,  which  produces 
as  output  a  standard  FORTRAN  progran.  which  any  compiler  can  then 


accept.  Alternatively, 


a  text  editor  nay  be  used  to  perforn  the  text 
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substitutions  inplied  t«v  the  Macros. 


In  the  work  that  is  described  here,  a  simple  Macro  preprocessor  based  on 
that  of  Kermghan  and  Plauger  (1976)  was  used. 


The  following  predefined  Macros  are  assuMed: 


*NACRO(NAME, VALUE)  Define  a  new  Macro. 

♦  INCLUDE ( f 1 lenane )  Include  a  file  in  the 

source  code  at  this 
point. 

tDELTOK  Delete  the  next  input 

token,  nornally  a  new 
line  character. 


The  reMaining  Macros  are  defined 
uses: 


SNAX 

5  i  ze 

IVHAX 

Size 

ICMAX 

Size 

SERAS 

Size 

in  the  code  and  have  the  following 

of  array  S. 
of  array  IV. 
of  array  IC. 
of  array  ERAS. 
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NAPMAXH 
HAf'MAX  V 
OUTLUN 
CHARACTER 


Horizontal  pii ge  size  in  character  widths. 
Vertical  page  size  in  character  heights. 
Logical  unit  nunber  for  the  output  device. 
Type  for  character  datim. 
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.5  T  HL  LjR  <  D 

The  computational  inJ  consists  of  an  irregular  triangular  tesselation 
which  is  tailored  to  approximately  correspond  to  the  shape  of  the  ocean 
basin  and  any  islands  that  are  to  be  modelled.  In  order  to  be  able  to 
conveniently  handle  ;uch  a  grid  on  a  computer,  we  must  formulate  a  data 
structure  that  represents  the  connect  .vitv  and  the  geometry  of  the  grio. 
First  let  us  introduce  two  definitions. 

A  vertex  is  defined  as  the  point  of  intersection  of  a  number  of  grid 
lines. 

A  cell  is  a  triangular  region  bounded  by  three  grid  lines. 

Next  we  must  rpfer  specifically  to  the  code  (given  as  an  appendix)  in 
order  to  See  row  the  data  structure  i  or g.imsed.  Cell  connectivity 
information  is  field  in  the  array  IC  and  vertex  connectivity  information 
is  held  in  the  array  IV. 

An  integer  expression  I  is  said  to  point  to  a  tell  if  IC<1)  is  the  first 
word  describing  the  cell,  similarly  an  expression  J  is  said  to  point  to 
a  vertex  if  IVt'J)  is  the  first  word  describing  the  vertex. 

In  addition,  a  number  of  words  of  real  storage  are  associated  with  each 
cell  and  verte:  for  storing  physical  quantities.  These  are  held  in  the 


t • 

ii 
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array  S. 


The  first  two  words  of  storage  for  vertices  ere  assumed  to  be  the 
coordinates  of  the  vertex.  Uhere  a  vector  field  is  to  be  represented, 
adjacent  words  in  S  are  used  to  hold  the  components  of  the  vector. 

As  an  example,  consider  the  shallow  water  equations.  The  physical 
storage  layout  associated  with  vertices  is  given  below: 

Uord 

0,1  Coordinates  of  the  vertex. 

2-8  Ueights  for  lumping. 

9-15  Transport  weights  in  the  x-di recti on. 

16-22  Transport  weights  in  the  y-direction. 

23  Depth. 

24,25  Momentum. 

26-23  Momentum  flux. 

29-31  Lumped  depth  and  Momentum  field  A. 

32-34  Lumped  depth  and  momentum  field  P. 

35  32  Lumped  depth  and  momentum  field  C. 

38.39  Force. 

40,41  Coriolis  term. 

The  meanings  of  these  fields  are  described  elsewhere. 
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The  layout  of  connectivity  information  for  a  cell  is  as  xuven  below: 

1 C ( I )  Pointer  to  the  S  array. 

lt'1+1)  Pointer  to  the  1st  vertex  that  defines  the  cell. 

1CU  +  2)  Pointer  to  the  2nd  ver((x  that  defines  th(  cell. 

IC(I+3)  Pointer  to  the  3rd  vertex  that  defines  the  cell. 


The  layout  of  connectivity  information  for  a  vertex  is  ■jiven  below: 


It'd  ) 
IVU  +  I ) 
I V  ( I  +  2  ) 
I V  <  I  >  3  ' 
I V ( I  +  4  ) 
I V  ( I  ♦  5 ) 
1VU+6) 
IV< 1+7) 


Pointer  to  the  S  array. 

Number  of  adjacent  vertices,  H. 

Pointer  to  the  1st  adjacent  vertex. 
Pointer  to  the  2nd  adjacent  vertex. 
Pointer  to  the  3rd  adjacent  vertex. 
Pointer  to  the  4th  adjacent  vertex. 
Pointer  to  the  5th  adjacent  vertex. 
Pointer  to  the  6th  adjacent  vertex,  or 
link  to  next  boundary  vertex. 


The  qrid  is  assumed  to  have  the  same  topoloqv  as  a  tessellation  of 
equilateral  triangles.  Arbitrary  boundaries  may  be  imposed  subject  to 
the  constraint  that  interior  boundaries  (islands!  do  not  have  acute 
anjles.  Since  the  qrid  is  deformed  before  use,  this  does  not  impose 
any  restriction  on  the  actual  shapes  of  islands  that  may  be  handled, 
however  the  qrid  would  be  more  than  usually  deformed  in  the  neighborhood 


J 


1U 


of  such  acuto  features. 

At  this  point  it  is  appropriate  to  examne  portions  of  the  code  in 
detail,  so  ue  will  examne  the  subroutines  that  plav  an  mportant  part 
in  constructinq  the  grid  data  structure. 

HEXTOP  is  a  subroutine  that  constructs  a  grid  with  the  sane  topology  as 
a  hexagonal  region  divided  into  triangular  cells.  It  is  called  by: 

CALL  HEXTQP(N) 

Each  side  of  the  hexagon  is  divided  into  N  segnents. 

f v  inspection,  we  nay  see  that  the  grid  has  6  N  triangular  cells, 
3(3N+1)N  edges  and  3(N+1)N+1  vertices. 

HEXTOP  generates  the  pointers  to  vertices  that  are  stored  in  IC,  by 
using  the  subroutines  GUC  and  GDC  which  generate  upward'  and  ■'downward' 
pointing  triangular  cells.  HEXTOP  counts  the  total  nunber  of  ceils 
allocated  in  NUMC,  and  the  nunber  of  vertices  in  NUMV. 

Next,  HEXTOP  generates  the  links  in  each  vertex  data  structure,  by 
examining  each  cell  and  Making  sure  that  all  vertices  surrounding  a  cell 
are  joined  together  by  using  the  subroutine  VLINK. 
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Finally,  HE x T OP  determines  which  vertices  lie  on  the  boundary  of  the 
region  tv  co-.intii.-i  the  number  of  adjacent  vertices.  The  boundary  points 
are  contained  in  a  linked  list  data  structure. 

HEXTOP  is  designed  as  a  specific  tool  for  generating  a  class  of 
topologies  that  are  useful  for  investigating  the  properties  of  advection 
schemes  and  in  simulations  of  anv  unbroken  region  that  has  approximately 
circular  boundaries.  In  addition  HEXTOP  illustrates  the  methods  that 
may  be  used  to  generate  even  more  complicated  topologies,  such  as  grids 
with  imbedded  holes,  from  the  human  engineering  standpoint,  an 
interactive  computer  system  with  a  display  screen  and  light  pen  would  be 
more  convenient  for  generating  grids  from  real  ocean  maps. 

GDC  and  GUC  are  cell  generating  subroutines  for  downward  and  upward 
pointing  cells  in  a  triangular  tesselation.  A  downward  pointing  cell  is 
generated  bv  the  call: 


CALL  GDC (IP) 


where  IP  points  to  the  upper  left  hand  vertex.  It  is  assumed  that  the 
following  vertex,  IP+IVSIZ.  is  the  upper  nqht  hand  vertex  and  that 
the  most  recently  defined  vertex  is  the  lower  verte:.  After  generating 
the  cell-vertex  links  for  a  new  cell,  GDC  returns  after  incrementing  IP 


to  the  next  vertex. 


GUC  works  rn  a  comp lementary  fashion. 


A  C  a  i  1 

CALL  GUC (IP) 

assumes  that  IP  points  to  the  upper  verte  ot  a  cell.  T Pie  nost  recent 
vertex  is  assorted  co  be  the  loner  left  vertex,  so  GUC  generates  a  neu 
verte;:  for  the  lower  right  hand  vertex  and  constructs  a  neu  set  of 
cel  1 -vertex  pointers.  IP  is  left  unchanged. 

By  alternating  calls  to  GUC  and  GDC  it  is  simple  to  construct  an 
arbitrary  tesselation  of  triangular  cells,  with  all  the  required 
pointers  to  define  the  topology  of  the  grid. 

VLINK  is  a  subroutine  that  ensures  that  two  vertices  are  linked  together 
in  the  order  given.  It  is  called  as  follows: 

CALL  VL I  NK  (  IVA,IVB) 

Verte-.  IVA  is  examned  to  see  if  it  is  already  narked  as  having  I VB  as  a 
neighbor,  if  so  VLINK  exits  without  doing  anything  extra.  Otherwise 
the  vertex  count  for  IVA  is  incremented  and  a  link  to  1VB  is  placed  in 
th.  next  available  vertex  link  p  sition.  1 o  correctly  link  two 


vertices,  two  calls  to  VLINK  are  required: 
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C-M.L  VLINKd'Jfl.IVPi 
CALL  VLINKUVB.IV'.' 


The  subroutine  YSTQ  assigns  storage  space  to  the  vertices  fron  the  array 
of  storage  S.  It  is  called  as  follows: 

CALL  VST0( IORIG,NUDS) 

The  first  word  allocated  is  = ( I  OR  16)  and  a  total  of  NUDS  words  are 
assigned  for  each  verte:;.  NUDS  should  be  equal  to  the  mmber  of 
physical  variable  fields  required  in  the  solutions  of  the  equations  of 
hydrodynamics. 

CIRBND  is  an  example  of  a  procedure  that  fixes  the  co-ordinates  of  the 
boundary  vertices,  in  this  case  by  distributing  them  uniformly  in  a 
circle.  It  may  be  called  without  any  arguments. 

CIRBND  is  used  with  HEXTOP  and  SUGRID  t  define  a  circular  ocean  basin 
for  testing  purposes. 

An  assumption  is  made  in  CIRBND  that  all  boundary  vertices  lie  on  a 
single  exterior  boundary,  and  that  the  links  between  boundary  vertices 


follow  serially  around  the  boundary.  This  is  true  for  HEXTOP,  but  need 
not  be  so  for  other  topology  generation  subroutines.  A  more  general 
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boundary  fixing  routine  would  be  needed  for  Multiple  boundaries,  and 
shouldprobabl v  be  enbedded  ir.  the  franework  of  an  interactive  computer 
svsten.  as  discussed  elsewhere. 

VCPUhR  is  a  subroutine  for  printing  out  the  links  associated  with  each 
vertex  and  cell.  It  is  called  without  anv  argunents.  This  routine  has 
pr-oved  useful  as  a  diagnostic  during  the  debugging  phase  of  new  topology 
generation  subrout i .its.  It  is  also  invoked  when  code  detects  an 
inconsistency  in  the  links,  caused  by  a  logic  error  m  the  topology 
routine,  or  More  usually  by  storage  corruption  caused  by  array  bound 
overflow  or  subroutine  argunent  inconsistencies,  which  traditionally 
are  not  detected  by  FORTRAN  conpilers. 

The  last  subroutine  in  the  grid  building  suite  is  SUGRID,  which 
calculates  the  co-ordinates  of  all  the  interior  vertices.  It  is  called 
as  follows: 


CALL  SUGRIIK NITER) 

Interior  vertices  are  Moved  so  that  the  co-ordinates  of  each  vertex  is 
equal  to  the  average  of  those  of  all  its  neighbors.  This  involves 
solving  a  svsten  of  2N  linear  equations,  where  N  is  the  number  of 
interior  vertices. 


The  equations  are  anenable  to  a  relaxation  process  of  the  wost 


straightforward  kind.  A  relaxation  parameter  REIIA  is  used  to  improve 
the  convergence  of  the  process. 


Empirical  tests  with  simple  grids  of  various  sizes  have  shown  that  RELPA 
=  1.388  seems  to  be  close  to  the  optimum  value. 

A  total  of  NITER  iterations  is  performed,  and  for  each  iteration,  the 
maximum  displacement  of  a  vertex  is  printed,  in  order  to  illustrate  the 
convergence  of  the  iterative  solution. 

On  the  basis  of  experience  with  other  relaxation  methods  applied  to 
elliptical  systems  of  partial  differential  equations,  it  is  believed 
that  SUGRID  is  unconditionally  stable  for  a  finite  range  of  RELPA. 

The  suite  of  subroutines  that  has  been  described  above  has  been 
specifically  written  to  construct  the  data  structure  for  a  circular 
ocean  basin,  but  simple  modification  to  CIRBNP,  for  example  would 
permit  irregular  quasi-c  n  cular  basins  to  be  treated.  If  the  basin  was 
grossly  dissimilar  to  a  circle,  or  had  a  different  topology,  for 
example,  containing  islands,  then  the  routine  HEXTDP  would  have  to  be 
replaced.  The  following  algorithm  is  proposed  for  the  generation  of 
grids  for  real -wo rid  oceanographic  simulations: 

First  assemble  a  number  of  equilateral  triangles  to  approximately 
represent  the  region  of  interest.  For  example,  HE  X  T  OP  uses  six 
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1  b 


triangles  in  the  torn  of  a  hexagon  to  represent  a  circle. 

The  next  step  is  to  subdivide  the  triangles  into  a  tesselation  of 
^nailer  equilateral  triangles  until  it  is  estimated  that  there  are 
sufficient  cells  to  resolve  the  solution  structure  that  is  desired. 

Triangular  cells  are  next  renoved  fron  the  interior  to  represent,  islands 
and  frow  the  exterior  to  better  represent  the  ocean  shore.  The  vertices 
on  the  exterior  and  interior  boundaries  are  then  assigned  the  actual 
coordinates  of  associated  points  on  the  actual  coastline. 

The  final  step  is  to  relax  the  interior  points  to  obtain  a  grid  with  a 
srtooth  transition  of  cell  sire.  The  whole  procedure  would  best,  be  done 
on  a  computer  systen  under  interactive  connand.  A  video  graphics 
display  and  a  lightpen  would  be  dost  appropriate  for  adjusting  the 


coastline  points  in  order  nate  the  grid  as  uniform  as  possible. 


4  IKE  ALGDftl  I  HU  I  UIV  NY CftODY AMI CS 

in  nost  fluid  dvnamcs  simulations  the  nost  troublesone  terns  in  the 
equations  are  the  non -linear  ones,  (he  advection  of  a  scalar  field  nay 
he  used  illustrate  the  algorithm  that  is  used  for  More  conpl icated 
cases. 


The  governing  equation  for  a  scalar  field  is: 


+  V  .F  -  0 

M 


Alternatively,  bv  Stoles  theoren,  we  can  say  that  for  any  region  1, 
hounded  hv  a  curve  5,  we  have 

i 

I  I 

♦  \  L-±-  z  ■' 

"\ 

where  ds  is  nornal  to  the  cm  ve  S,  and  ^  is  the  average  value  of  ^ 
in  the  region  1 . 

In  a  finite  representation  of  <f>  using  a  grid  of  points  on  a  triangular 
nesh,  it  is  convenient  to  store  the  value  of  y>  at  the  vertices  of  the 
grid,  and  to  consider  a  region  surrounding  each  vertex  which  we  will 
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coll  a  flux  cell. 

A  flu.,  cell  is  on  irregular  and  sonetines  non-conve>:  polygon  that 
contains  only  one  vertex.  Each  side  of  the  polygon  starts  at  the 
cent"oid  of  a  cell  adiacent  to  the  contained  vertex  and  finishes  on  the 
md-point  of  a  side  of  the  cell  that  passes  through  the  contained 
vertex.  It  is  clear  that  no  part  of  the  region  defined  by  the  grid  lies 
outside  of  all  flux  cells. 

A  flux  cell  surrounding  a  vertex  with  H  neighboring  vertices  is  a 
polygon  with  2N  sides. 

Averages  of  a  variable  within  a  flux  cell  are  tented  limped'  values  and 
are  associated  with  the  contained  vertex.  If  the  values  of  the  variable 
are  to  be  represented  by  its  values  at  the  vertices  only,  as  is  usual 
for  anv  finite  difference  schene,  then  we  ftav  uniquely  approxinate  the 
variable  at  any  point  by  linear  interpolation  between  the  vertices  of 
the  cell  containing  the  point..  Uith  this  representation  of  the 
variable.  we  way  exactly  calculate  limped  values  bv  integrating  over 
the  flux  cell.  Similarly,  the  line  integral  above,  nay  be  exactly 
evaluated  around  the  flux  cell. 

Our  advection  equation  is  conservative  in  the  limped  values  pf  the 
scalar  variable.  Fluxes  are  derived  bv  first  evaluating  the  unlimped 
values  at  vertices,  evaluating  the  non-linear  flux  terns  at  each  vertex 


_ 


and  then  rer forming  linear  interpolation. 


Thus,  we  have  a  simple,  non-anbiguous  algorithm  for  discretising  the 
continuous  equations  of  hydrodynamics. 

The  line  and  surface  integrals  of  interpolated  fields  nay  be  exactly 
representeo  by  a  weighted  sun  of  vertex  values  in  the  vicinity  of  the 
region  of  integration.  The  weights  do  not  change  with  tine,  so  the 
complicated  integrals  may  be  simply  evaluated  with  a  minimum  of 
computational  effort. 

The  subroutine  CVLU  calculates  the  weights  A,  Bj  required  tD  calculate 
the  vertex  lumped  values  Fi 

Fi  =  A  Fi  ♦  £  F.i  Bi 

where  ,i  is  a  vertex  which  neighbors  vertex  1. 

Consider  a  vertex  i  and  a  cell  k.  The  field  F  1,  defined  at  l  and  the 
other  two  vertices  p  and  q  such  that  the  cell  f  is  defined  by  the 
triangle  lpq.  If  we  approximate  F  within  F  by  linear  interpolation,  we 
may  integrate  F  within  the  quadrilateral  ip'gq'  where  q  is  the  midpoint 
of  ip  and  p  is  the  midpoint  of  ip.  and  g  is  the  centroid  of  ipq. 

This  integral  is 


■IWdOOBi 


V*  =.»k  I  11/54  1-1  +  7/108  Fp  +7/108  F.qJ 


where  K  is  the  area  of  cell  k. 

The  vertex  limped  value  of  F  at  l  is  the  sun  of  over  all 
quadrilaterals  as  above  which  have  l  as  a  vertex.  The  weights  A  and  Bj 
flay  thus  be  seen  to  be 


A  +  -11/54  y  \  k. 

x. 

V  . 

B.i  =  7/108  '  A  Ik 

t- 

where-^jk  3  4  k  if  j  is  a  vertex  of  k,  else  L.  jk  -  0. 

CVLU  is  invoked  as  follows: 

CALL  CVLUl IPOS) 

IPOS  is  a  pointer  to  a  set  of  7  variable  locations  which  are  to  contain 
the  values  of  A  and  Bj. 

The  subroutine  CTU  conputes  the  weiqhts_A  and  Pj  to  calculate  the  line 
integral  T  of  a  vector  flux  F  around  a  flux  cell  surroundinq  a  vertex  l 
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where 


and  .1  is  a  vertex  neighboring  1. 

A  and  B j  are  vectors  associated  with  vertex  1  and  are  stored  with  their 
x-conponents  at  IPOS  and  IPOS+J  respectively.  The  y-conponents  are 
stored  7  locations  after  the  corresponding  x-conponent. 

The  flux  cell  around  a  boundary  vertex  is  truncated  by  the  actual 
boundary,  so  that  part  of  the  line  integral  Must  be  evaluated  along 
sections  of  the  boundary,  however,  straightforward  use  of  linear 
interpolation  suffices  to  calculate  A  and  Bj.  The  details  of  the 
algorithn  need  not  be  described  here,  since  the  code  itself  serves  to 
define  the  Method. 

The  subroutine  is  called  as  shown: 

CALL  CTU(IPOS) 

where,  as  before,  IPOS  is  a  pointer  to  the  region  of  storage  that  is 
to  contain  the  weights.  In  this  case  14  words  are  required  to  store  the 


1 


vector  weights. 
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The  correctness  of  the  algorithm  nay  he  tested  by  the  subroutine  TSTCTU, 
which  evaluates  the  line  integrals  for  a  flux  of  (1,0)  and  (0,1)  for  the 
boundaries  of  the  flux  cells  containing  each  vertex.  Within  the  linits 
of  rounding  errors,  the  two  integrals  are  zero. 

TSTCTU  is  called  as  shown: 

CALL  TSTCTU(IPOS) 

where  IPOS  is  the  pointer  to  the  weights  generated  by  CHI. 

LUMP  is  the  subroutine  that  calculates  the  integral  of  a  physical 
variable  field  over  the  flux  ceil  that  surrounds  each  vertex. 

The  subroutine  LUMP  is  invoked  as  follows: 

CALL  LUMP ( IPOS , I  FROM , I  TO  ) 

IPOS  is  a  pointer  to  the  set  of  weights  that  have  previously  been 
calculated  by  CVLU,  and  IFROM  and  ITO  are  pointers  to  the  source  and 
destination  fields  respectively. 

The  limped  variable  is  stored  at  location  ITO  after  calculating 
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where  A  and  Bi  are  the  weights  and  Fj  is  the  value  of  the  field  at 
vertex  j.  j  is  a  neighbor  vertex  of  vertex  1. 

UNLUMF  is  a  subroutine  which  is  the  formal  inverse  of  LUMP'.  It  will 
recover  the  original  values  of  a  field  that  has  previously  been  smoothed 
by  LUMP.  It  is  invoked  as  shown: 

CALL  UNLUMP ( IPOS, I FROM, I  TO, ERR ) 

IPOS,  IFROM  and  ITO  have  the  sane  meaning  as  in  the  above  description 
of  the  subroutine  LUMP. 

The  systen  of  linear  equations  for  Fi,  given  Fi  is  solved  by  relaxation 
Methods,  with  a  final  accuracy  estimated  to  be  on  the  order  of  ERR  if 
possible,  otherwise  an  error  Message  is  printed. 


1  he  subroutine  ABVECT  is  called  as  follows: 

CALL  ADVECT  < IPCN.IPCL, IPC, IPU, IPU.DT  ) 

This  subroutine  perforMs  the  basic  advection  of  a  limped  scalar  field  at 
location  IP CL  using  the  unlunped  velocity  whose  x-  and  v-coMponents  are 
located  at  IPU  and  IPUH,  Transport  weights  are  located  at  IPU.  The 


tine  interval  for  integration  is  DT  and  the  new  scalar  field  is  returned 


to  the  location  IPCN. 


The  subroutine  ADVECT  Modifies  a  lumped  rtonentun  field  to  ensure  that 
the  corresponding  momentum  values  on  the  grid  boundaries  are  zero.  This 
is  equivalent  to  imposing  a  rigid  boundary  condition  at  the  edge  of  the 
computational  domain. 

The  subroutine  STEP  is  called  as  follows: 

CALL  5TEPUPLN ,  IF'ULN ,  II-'VLN ,  IF'LO,  IPULO,  IPVLO, IPL, 
IPUL. 1PVL. IP, IPU, IPV,IPUU,IPUV,IPVV, 
IPFX,IPFY,IPFXL,IPFYl,IPU,DT,F,CF, 

USX , USY ) 

The  STEP  subroutine  advances  the  fields  of  the  physical  variables  one 
time  step  for  the  shallow  water  equations.  The  values  of  the  dependent 
variables,  height  and  horizontal  momentum  components,  at  the  start  of 
the  step  are  denoted  with  a  suffix  '0'  for  Did  while  the  values 
corresponding  to  the  end  of  the  step  have  a  suffix  '  N'  for  '’New'.  STEP 
first  unlumps  the  fields  in  order  to  calculate  pressure  and  Reynolds 
stress  terms,  then  uses  the  transport  algorithm  and  finally  imposes  the 
boundary  conditions. 


The  subroutine  TRNSPT,  which  is  called  as  shown: 
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CALL  TRNSP'f  <  INEU,  IOLU,  1FX,  IFY,  IPU.I'T  ) 

evaluates  the  changes  to  a  lunped  field  for  a  tine  step  DT.  The 
original  lunped  field  is  located  at  IOLD.  The  field  after  Modification 
for  the  effects  of  transportation  is  returned  to  location  INEU.  The 
transportation  is  effected  by  a  flux  field  whose  x-  and  y-conponents  are 
located  at  IPX  and  IFY  respectively. 

The  shallow  water  equations  contain  several  terns  which  are  not  of  an 
advective  nature.  These  terns  are  generally  easier  to  treat  in  a 
numerical  schene  than  the  non-linear  advective  terns.  Ue  will  refer  to 
the  non-advective  terns  as  forcing  terns.  The  subroutine  FRF  includes 
the  effects  of  the  following  forcing  terns: 

1)  Frictional  tern  with  a  coefficient  of  friction 
Cf  and  a  linear  dependency  with  the  fluid 
speed. 

2)  Coriolis  tern  with  a  paraneter  F. 

3)  An  externally  inposed  stress  (USX.USY)  which 
nay  be  regarded  as  representing  the  effects  of 
wind  on  the  fluid. 


f 


I 

\ 


The  subroutine  u  called  with  the  arguments 


FRF  ( IF'L  f  IPUL ,  IPVL  ,  IP,  IPU.IPV,  IF'FX  ,  IPFY  ,  IF'FXL  , 
!FFTL.DT,F,CF,USX,«Sr) 

IF',  IF'U  and  I F V  ore  pointers  to  unlumped  depth  and  momentum  components. 
I  PI. ,  IPUL  and  I P  V  L  are  the  locations  of  the  corresponding  lumped 
fields.  IF'FX  and  IF'FY  die  fields  used  to  store  the  combined  frictional 
and  external  stresses.  IFF  XL  and  IPFYL  are  used  to  locate  the  lumped 
values  of  the  previous  two  fields.  Upon  exit  from  FRF  the  momentum 
fields  are  modified  to  reflect  the  time  integrated  effects  of  the 
forcing  terms  over  a  time  interval  DT, 

Two  main  program  are  included  in  the  appendix.  The  first  is  used  in 
simple  test  solutions  for  the  color  equation,  in  order  to  examine  the 
stability  and  accuracy  of  the  basic  advection  algorithm. 

The  second  subroutine  is  a  driver  for  the  shallow  water  equations  and 
includes  calls  to  the  graphics  routines  which  are  described  in  the  next 
section.  The  subprogram  makes  two  calls  to  STEP,  since  the  basic  time 
differencing  scheme  that  is  employed  is  the  robust  pseudo-backward  Euler 


method . 
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'j  UU1PUT  ROUTINES 

In  debugging  the  code  and  evaluating  the  perf ormanc e  of  the  numerical 
scheme  that  it  embodies,  graphical  output  is  an  invaluable  tool. 
Simple  listings  of  field  values  give  little  indication  of  the  actual 
behavior  of  the  Model  being  simulated,  especially  when  an  irregular 
grid  is  used.  A  number  of  routines  have  been  produced,  which  serve  m 
themselves  as  valuable  aids,  while  at  the  sane  tine  forning  a  set  of 
pnmtives  for  nore  elaborate  displays. 

The  QPRCON  subroutine  is  a  sinple  display  procedure  for  producing 
contour  plots  of  a  field  on  a  hardware  device  such  as  a  lineprmter.  It 
is  called  as  follows: 


CALL  QPRCON!IFLD,XNIN,XNAX> 


I FLD  is  a  pointer  to  the  field  that  is  to  be  displayed.  A  single 
character,  either  a  digit  or  a  plus  or  Minus  sign  is  placed  on  the  page 
at  a  position  corresponding  to  a  vertex  location.  The  character 
represents  the  value  of  the  specified  field,  with  0  corresponding  to 
XN  IN  and  9  corresponding  to  XNAX.  Intermediate  values  are  represented 
by  an  appropriate  decinal  value  using  a  linear  transfer  function. 
Values  less  that  XNIN  are  shown  by  a  Minus  sign  while  values  larger  that 


XNAX 

are  shown 

with 

a  plus  sign. 

If  the  grid  is 

coarse,  the  display 

will 

be  sparse, 

but 

QPRCON  has  the 

advantage  of 

representing  both  the 

grid  structure  and  the  field,  so  that  lodgements  mav  be  made  as  to  the 
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reality  of  features  with  scales  near  that  of  the  basic  grid. 

PRCON  is  a  subroutine  that  is  functionally  sinilar  to  QPRCON,  having 
the  sane  argument  list,  except  that  the  field  values  for  space  between 
vertices  are  also  displayed.  Each  character  space  on  the  output  page  is 
napped  into  the  grid,  and  if  the  point  lies  within  the  boundary  of  the 
grid,  the  value  of  the  input  field,  IFLD,  is  evaluated  by  linear 
interpolation.  The  running  tine  for  PRCON  nay  be  nurh  larger  than  that 
for  QPRCON,  especially  when  the  the  nunber  of  vertices  is  less  that  the 
nunber  of  resolvable  positions  on  the  output  page. 

The  routine  CONPLT  generates  a  contour  line  plot  of  the  specified  field. 
It  is  called  by  neans  of: 

CALL  CONPLT (IFLD,C ON VAL,N) 

Contour  values  are  printed  for  values  CONVAL(I),  1=1 ,N.  CONPLT  is 
suitable  for  all  types  of  contour  plotting  devices  such  as  pen  plotters 
and  electrostatic  printer-plotters,  since  it  nak.es  use  of  only  one 
svsten  dependent  subroutine,  LINE,  for  drawing  a  straight  line 
segnent . 

The  DRAUV  subroutine  is  called  as  follows: 


DRAUV(IPL,IPUL, IPVL , IP, IPU, IPV , ASCALE, I PU ) 
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The  DRAUV  subroutine  generates  a  display  of  the  velocity  field  implied 
bv  the  lumped  depth  and  momentum  fields  located  at  IPL,  IPUL  and  IPVL. 
The  fields  IP,  IPU  and  1PV  are  used  as  workspace  for  calculating 
unlumped  values  of  the  physical  variables. 

The  COMMON  variable  SCALE  is  used  to  transform  from  grid  coordinates  to 
plotter  coordinates.  After  unlumping  the  fields,  the  velocity 
components  are  extracted  from  the  momentum  fields.  At  each  vertex  of 
the  grid,  an  arrow  is  drawn  with  a  length  proportional  (using  the  scale 
factor  ASCALE)  to  the  velocity  at  the  specified  point.  No  attempt  is 
made  to  draw  arrows  that  are  so  small  as  to  be  unresol vable. 

BRGRIP  is  a  subroutine  that  has  no  arguments  but  that  may  be  used  for 
drawing  the  triangular  grid  on  a  scale  similar  to  the  other  plotting 
displays,  so  that  it  may  be  used  to  overlay  contour  dt  vector  field 
representations. 

The  subroutine  DRCELL  is  similar  to  DRGRID,  also  having  no  arguments, 
e.cept  that  it  draws  the  outline  of  the  flux  cells  that  are  used  for  the 
line  integrals  in  Stokes'  theorem. 

EUPLOT  is  a  subroutine  with  the  following  calling  sequence: 


CALL  EUPLOT(IP.ZMAX,ZMIN) 
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The  EUPLOT  subroutine  is  used  to  generate  a  clot  of  the  field  located  at 
IP,  along  the  line  T  =  0.  The  actual  plot  is  produced  by  printing 
characters,  in  a  way  suitable  for  use  with  a  linecrinter.  The  «a;:iMUM 
and  minimum  values  that  nay  be  displayed  on  the  page  are  given  by  ZNAX 
and  ZMIN  respectively. 

INTPOL  is  the  basic  interpolation  subroutine  used  by  the  other  graphics 
subroutines  for  extracting  field  values  fron  arbitrary  positions  on  the 
grid  bv  means  of  linear  interpolation.  It  is  called  as  follows: 

CALL  INTPOL (IFLO,X,Y,Z,FAIL) 

The  postion  specified  is  given  by  the  coordinates  (X,T)  and  the  value  of 
the  field  indicated  by  IFL  is  returned  in  the  variable  Z.  If  the  point 
lies  outside  the  boundary  of  the  grid,  the  logical  variable  FAIL  is  set 
to  betrue.  INTPOL  works  by  examining  each  cell  in  turn  until  one  is 
found  that  contains  the  required  point,  so  it  is  relatively  tine 
consuming,  especially  when  it  is  called  manv  tines,  as  in  the 
subroutine  PRCON. 

The  routine  could  be  increased  in  speed  by  introducing  a  relatively 
coarse  rectangular  grid  and  associatinq  a  linked  list  of  triangular  cell 
numbers  with  each  rectangular  cell.  The  search  could  then  be  performed 


on  a  much  shorter  list,  the  particular  list  being  determined  by  direct 


r 
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computation  from  X  and  1.  This  would  he  analogous  to  hash  coded 
searches  in  standard  table  look  up  problems,  with  the  added  advantage 
that  each  list  would  be  approximately  the  same  size,  so  the  look  up 
time  would  be  both  short  and  predictable. 

INTRI  is  function  that  return  a  value  of  LOGICAL  type.  It  may  be 
referenced  by: 


INTRI ( X , Y ,  ICELL  ) 

The  returned  value  is  only  true  if  the  point  (X,Y)  lies  inside,  or  on 
the  boundary  of  the  triangular  cell  number  ICELL. 

The  subroutine  MAXHIN ,  which  is  called  by: 

CALL  MAXNINUFLD) 

prints  out  the  maximum  and  minimum  values  of  the  field  specified  by 
IFLD. 


The  function  NEIGH  returns  a  LOGICAL  value  if  the  two  vertices  specified 
by  IP  and  10  are  neighbors  when  it  is  called  by  the  sequence: 


NEIGH( IP, IQ j 


32 


The  subroutine  PRFLDS ,  called  by 

CALL  PRFLDSIII ,12) 

prints  out  the  values  of  all  the  fields  in  the  range  fron  II  to  12  in 
tabular  forn  for  each  vertex.  If  needed  the  table  is  split  into  several 
parts  so  as  not  to  exceed  the  width  of  the  lineprinter  page. 

The  VCDUMP  subroutine  displays.  on  the  printer,  all  the  links 
associated  with  the  cell  and  vertex  data  structures.  VCDUHP  does  not 
take  any  arguments. 

CONSRV  is  a  subroutine  called  as  follows: 

CALL  CONSRVT IFROW) 

This  diagnostic  routine  evaluates  and  prints  the  sun  of  the  values  of  a 
variable  at  each  vertex  of  the  grid.  If  the  variable  is  a  lumped  field, 
the  sun  is  the  spatial  integral  of  the  corresponding  unlunped  field. 
The  routine  has  proved  itself  useful  in  code  testing  by  evaluating  the 
integral  of  fields  that  should  be  fornally  conserved  in  tine  by  the 
systen  of  equations. 

The  COPY  utility,  which  is  called  by: 
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CALL  COPY ( I  TO, IFROM) 

is  used  in  nany  places  throughout  the  code  to  transfer  a  physical 
variable  field  fron  one  location  IFROM  to  another  location  ITO. 

The  subroutine  FLDAV  is  a  utility  procedure  uhich  is  called  by: 

CALL  FLDAV< IPOSA, IPOSB,  IPOSC) 

It  calculates  the  average  of  two  fields  indicated  by  IPOSB  and  IPOSC  and 
returns  the  result  to  the  location  IPOSA. 

INITF  is  a  subroutine  used  to  set  up  initial  fields  of  unit  depth  and 
zero  nonentun.  It  is  called  by: 

CALL  INI TF ( IP , IPU , IPV ) 

The  subroutine  5BR0T  inposes  a  velocity  field  appropriate  for  a  solid 
body  rotation  with  an  angular  velocity  of  OMEGA.  The  resulting 
components  are  stored  at  locations  IPOSU  and  IPOSU+1  when  it  is  called 
by : 


CALL  SBROKQMEGA, IPOSU) 


The  TFIEID  subroutine  sets  up  a  scalar  test  field  which  is  a  compact, 
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continuous  and  differentiable  piecewise  bivariate  cubic,  in  the  shape 
of  an  off-center  'hump'.  This  field  has  proved  useful  for  studying  the 
advection  of  a  scalar  field  under  solid  body  rotation.  Although  this 
problem  is  analytically  trivial  since  the  field  nerely  rotates  with  the 
flow  without  changing  shape,  it  is  quite  non-tnvial  computationally 
and  provides  a  powerful  test  for  comparing  the  efficacy  of  different 
numerical  advection  schemes.  It  may  be  invoked  by  the  call: 


CALL  TFIELD(IPOS) 
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6  RESULTS 

A  sequence  of  numerical  experiments  uere  performed  with  the  code  during 
its  development  in  order  to  provide  answers  to  the  following  questions: 

1)  Can  a  simple  algorithm  be  developed  which  is 
capable  of  solving  the  equations  of  hydrodynamics 
with  the  physical  variables  approximated  by  values 
on  an  irregular  grid? 

2)  Uhat  are  the  limits  of  the  stability  of  the 
numerical  scheme? 

3)  Uhat  is  the  accuracy  of  the  scheme,  especially  in 
the  treatment  of  the  dominant  non-linear  advective 
terms? 

4)  Uhat  computing  resources  are  required  to  use  the 
algorithm? 

5)  How  does  the  algorithm  compare  to  orthodox  finite 
difference  schemes7 

The  answer  to  question  1)  is  clearly  yes.  The  scheme  that  has  been 
described  is  certainly  simple  and  elegant.  Some  slight  complications 


are  inherent  in  using  an  irregular  grid,  but  all  the  geometrical  and 


topological  factors  are  resolved  once  at  the  start  of  a  simulation  and 
are  incorporated  in  the  weight  terns  and  the  cell  and  vertex  links. 
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During  execution  the  scheme  is  fully  explicit,  except  for  the  simple 
relaxation  scheme  that  is  used  for  'unlumping'  the  physical  variable 
fields . 

One  especially  satisfactory  feature  of  the  algorithm  is  that,  for  the 
treatment  of  spatial  derivatives  at  least,  there  are  no  arbitrary 
decisions  to  be  made  in  setting  up  the  difference  equations,  Host 
finite  difference  schemes  are  overdetermined,  in  the  sense  that  several 
different  representations  have  similar  spatial  order  of  accuracy,  with 
often  no  obvious  means  of  resolution.  Also  the  present  algorithm  may  be 
readily  generalised  to  more  complicated  situations. 

The  answer  to  question  2)  is  not  capable  of  analytical  solution,  so  we 
must  perform  experiments  and  make  cautious  inferences.  Host  finite 
difference  schemes  that  are  explicit  in  time  differencing  must  have  some 
constraint  on  the  maximum  time  step  in  order  to  satisfy  stability 
criteria.  If  these  stability  criteria  are  not  satisfied  then  a 
computational  mode'  solution  may  be  obtained.  These  frequently  exhibit 
high  spatial  frequency  oscillations  and  bear  no  resemblance  to  the 
solutions  of  the  underlying  differential  equations. 


Some  numerical  schemes,  such  as  the  Dufort-Frankel  scheme,  while  being 
explicit  in  time,  have  no  formal  constraint  on  the  time  step.  However 


they  achieve  stability  at  the  expense  of  introducing  increasing  amounts 
of  artificial  diffusivity  as  the  tine  step  is  nade  larger.  In  addition, 
even  schenes  that  nake  use  of  mplicit  tine  differencing,  such  as  ADI 
nethods  are  still  in  practice  lmited  by  the  need  to  naintain  adequate 
accuracy  in  the  tine  integration. 

A  umber  of  nunerical  tests  have  been  performed  both  for  the  scalar 
advection  code  and  the  shallow  water  equations.  For  sem -regular  grids 
the  naxinun  tine  step  for  stability  was  observed  to  be  inversely 
proportional  to  the  nean  inter-vertex  spacing,  while  for  nore  irregular 
grids  the  naxinun  stable  tine  step  was  found  to  be  nore  closely  related 
to  the  inverse  of  the  mninun  inter-vertex  spacing.  The  latter  result 
is  hardly  suprising  but  it  does  nean  that  when  constructing  an  irregular 
grid  to  follow  a  coastline  of  creat  conplexity,  it  pays  not  to  use 
gratuitously  snail  cells.  However  if  the  boundary  points  are  reasonably 
distributed  along  the  coast,  and  the  grid  topology  is  appropriate  to 
the  conputational  c'onain  the  use  of  SUGRIP  or  its  generalisation  will 
ensure  that  no  excessively  snail  inter-vertex  spacings  are  generated. 

Overall,  it  has  been  found  that  the  stability  behaviour  of  the 
triangular  grid  schene  is  qualitatively  and  quantitatively  sirtilar  to 
nore  conventional  rectangular  grid  finite  difference  schenes. 

Sone  answers  to  question  3)  have  been  obtained  by  solving  the  color 
equation  for  a  solid  body  velocity  field.  Such  an  equation  has  a 
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trivial  analytics!  solution  so  we  nay  investigate  the  accuracy  of  the 
schene  in  its  creatnent  of  the  all  important  advection  tern.  For  a 
circular  basin  of  unit  radius,  a  solid  body  rotation  law  and  an  initial 
field  deternined  bv  the  subroutine  IF1ELD,  good  results  were  obtained 
even  with  a  very  coarse  grid.  For  example  with  a  typical  inter-vertex 
spacing  of  .125,  the  bell  shaped  test  field  only  suffered  about  103 
error  after  a  conplete  circuit  of  the  grid.  Since  positivity  is  not 
ensured  bv  the  schene,  sone  snail  negative  and  positive  ripples  were 
observed,  especially  in  the  wake  of  the  “hunp",  but  these  always 
renamed  snail  in  nagnitude.  The  algorithm  has  second  order  spatial 
accuracy  so  in  the  linit  of  fine  grids,  error  estinates  nay  be  easily 
calculated. 

Question  4  nay  be  answered  both  on  theoretical  estinates  and  on  the 
results  of  the  nunerical  expennents  that  have  been  performed.  If  we 
take  as  a  neasure  of  the  grid  conplexity  sone  nunber  N,  for  exanple 
NUNOIV  as  used  in  the  code,  or  the  reciprocal  of  the  inter-vertex 
spacing.  Then  the  nunber  of  grid  points  for  which  calculations  must  be 
perforned  is  N**2.  The  solution  of  the  systen  of  linear  equations  for 
the  unlunpini  operation  requires  on  the  order  of  N**3  operations,  since 
it  is  iterative,  and  infornation  nust  be  propagated  over  the  whole  grid 
in  order  to  effect  a  solution.  For  each  tine  step,  the  computational 
effort  has  two  factors,  one  proportional  to  the  square  of  N  and  the 
other  proportional  to  the  cube  of  N.  For  noderate  values  of  N  up  to 


about  20, 


the  first  tern  is  doninant 


For  fine  scale  sinulations, 


however,  the  effort  of  the  UNLUMP  operation  will  dominate  the  solution. 
The  maximum  pernitted  tiMe  step  vanes  inversely  with  N,  so  a  conplete 
sinulation  will  take  a  tiMe  that  is  proportional  to  the  cube  of  N  for 
coarse  grids  and  the  fourth  power  for  sufficiently  fine  grids. 

There  is  no  clear  cut  answer  to  the  final  question.  Sone  cartesian  grid 
Methods  Maintain  N**3  tiMing,  and  so  for  very  fine  grids,  will  be  mo re 
efficient  than  the  irregular  triangular  grid  Method.  On  the  other  hand, 
irregular  boundaries  May  require  excessively  fine  grids  for  the 
cartesian  grid  Methods,  so  their  theoretical  edge  May  not  be  Maintained 
in  practice.  It  has  been  denonstrated  that  irregular  triangular  grid 
Methods  can  be  applied  to  the  equations  of  ocean  flow,  without  undue 
penalties  in  conputer  resources,  or  conplexity  of  the  final  code. 
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Appendix 


Typicsl  set  of  eacro  definitions 

*MACR0-;*ft*DELT0K]>* 
»MACR0(SMAX,6510r 
♦MACRO!  I VMAX  ,1736)"' 

$ MACRO! I  CM AX, 1668)  " 
♦MACR0(MAPHAXHf72)“ 

*  MACRO  (MAPMAXV.22)'' 
♦MACRO!SERAS, 1 ) * 
*HACR0(0UTLUNf5)- 
•MACRO! INLUN , 5) * 

♦  MACRO (CHAR AC  TER  .LOGICAL  *ir 


A2 


Mam  program  for  driving  the  scalar  advection  code 


$INCLUDE(CHAC.MAC)*D£LTOK 
TYPE  1 

1  FORMAT (  '  N,  DT' ) 

ACCEPT  *,N,DT 
URITE<0UTLUN,2)  Nf0T 

2  FORMAT ( '  N,DT',I4,F12.5) 

CALL  CTEST (NfDT) 

STOP 

END 

SUBROUTINE  CTEST (NUMDIV,DT) 

URITE<  t ,96)  NUMDI V, DT 
96  FORMAT  (  '  N ,  DT  ' ,  16, FH  .6) 

0MEGA=6. 283195 
NSTEPS- 1 . /DT 
CALL  HEXTOP(NUMDIV) 

CALL  VST 0 ( 1,30) 

CALL  CIRBNB 
CALL  SUGRID(20) 

IPLW=2 

CALL  CVLU(IPLU) 

IPU=9 

CALL  CTU(IPU) 

C  CALL  TSTCTU(IPU) 

CALL  GRDLEN 
IP  =  23 
IU=24 
I  V=25 
IPL=26 
IP N*27 
IUL=28 
I Vl=2? 

CALL  TFIELD(IP) 

CALL  SBROT ( OMEGA, III) 

CALL  LUMP  ( IPIM ,  IP  ,  I F’L  ) 

C  Remove  the  Cs  fron  the  next  three  lines  to  inpose  rigid  b.c. 
C  CALL  LUMP ( IPIU, 1U, IUL  ) 

C  CALL  LUMP ( I PIU, I V, I VL ) 

C  CALL  RIGBND<IUL,IVL,IU, IV) 

DO  10  1=1 , NSTEPS 

CALL  ADVECT (IPN,IPL,IP,IU,IPU,DT) 

CALL  COPY ( IPL  f IPN) 

CALL  UNIUHPUPLU,  IPL,  IP,  1  .E-5) 

IF  ( MOO ( 1,10)  .EQ.  0)  CALL  HAXMINUP) 

10  CONTINUE 

CALL  HAXNINIIP) 

RETURN 

END 
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SUBROUTINE  GRDLEN 
♦INCLUDE ( TRG .  DCL  ) ' 

J=1 

NSEG=0 
XNIN*1 .  E 1 0 
XMAX=0. 

XNEAN=0. 

DO  1  I=1,NUNV 
JP= J+1 
NV=IV( JP) 

JA= I V  < J) 

X=S( JA) 

r=s ( ja+ i ) 

DO  2  K=1 , NV 
JK= JP+K 
JC= I V ( JK ) 

JB=I V  < JC) 

XX=S( JB) 

YY=S  <  JB+1 ) 

D=SQRT<  <X-XX)**2MY-YY)**2) 

NSEG=NSEG+1 

XHEAN=XNEAN+D 

IF  (D  .GT.  XMAX  >  X«AX=D 

IF  (D  .LT.  XNIN)  XNIN«D 

2  CONTINUE 

1  J= J+IVSIZ 

XNEAN=XNEAN/NSEG 

3  FORMAT ( '  HIN,HEAN,HAX' , 3F12.5) 

UR ITE< 1 ,3)  XNIN, XNEAN, XMAX 
RETURN 

END 

SUBROUTINE  TFIELD(IPOS) 

♦  INCLUDE ( TRG, DCL  )  ‘ 

YC*0 . 

XC  =  1. 

RC=.? 

1  =  1 

DO  1  11=1 , NUhV 
K=IV( I ) 

X  =  S  ( K ) 
r=S(K+i ) 

R=SGRT ( <  X-XC)**2+ <  Y-YC )**2)/RC 
C=0. 

IF  <R  .LT.  1.)  C=(2.*R-3.)*R«R*1. 
S  <  K ♦ IPOS) *C 
1  I=I+IVSIZ 
RETURN 
END 


i 

I 
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Nam  proqran  for  the  shallot*  water  equations  code 


♦  INCLUDE  <  TMAC  .MAC  )*I»ELT0K 

DIMENSION  ERAS ( SERAS ) 
DIMENSION  CONVAL(IO) 
LOGICAL  PLOT 
PL0T=. FALSE. 

ASCALE=1. 

DT=.04 

CALL  INPUT ( 'DT _ ',0T) 

X  =  2 

CALL  INPUT ( 'NUMDI V'  , X ) 

NUMDIV=X 

X=1 01 

CALL  INPUT ( 'NSTEPS' ,X) 
NSTEPS=X 
X  =  1 0 

CALL  INPUT (  'NCONT  ' ,X) 
NC0NT=X 
DO  2  1=1,10 
2  C0NVAL(I)=I*.2-.1 
CALL  HEXTOP(NUMDIV) 

CALL  VSTOt 1,42) 

CALL  CIRBND 
CALL  SUGRIDI20) 

IPLM=2 

CALL  CULU(IPLU) 

IPU=V 

CALL  CTU(IPU) 

CALL  TSTCTU(IPU) 

IP  =  23 

IPU=24 

IPV=25 

IPUU=26 

IPUV=27 

IPVM=28 

IPLA=2? 

IPULA=30 

IPVLA=31 

IPLB=32 

IPULB=33 

IPVLB=34 

IPLC=35 

IPULC=3A 

IPVLC=37 

IPFX=38 

IPFT=3? 

IPFXL= AO 
IPFYL=A 1 

CALL  INITF<IP,IPU,IPV) 
CALL  LUMP! 2 , IP, IPLA) 
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CALL  LU«P(2,IPU, IPULA) 

CALL  LUMP ( 2, IPV, IPVLA) 

CF=0. 

C  F  =  1  . 

CALL  INPUT  DCF - ,CF) 

F  -0 . 

CALL  INPUTT  F .  ,F> 

USX-  .3 

usy  =  o. 

CALL  INPUT*  USX...  .USX) 

CALL  INPUT <  USY...  , US V ) 

IF  (PLOT)  CALL  PLOTS (ERAS ,900,5.5) 

DO  10  1=1 ,NSTEPS 

IF  ( HOD  (1-1  .NCOfll  )  .NE.  0)  GU  TO  9 
CALL  UNLUNP<2, 1  PL A. IP, I .E-5) 

CALL  QPRCON< IP,0. ,2. ) 

CALL  EUPLOT ( IP,2. ,0. ) 

IF  (PLOT?  CALL  0RIGIN(5.5,0. ) 

IF  (PLOD  CALL  DRAUV ( IPLA, I PUL A , IPVLA, IP, IPU, IPV, ASCALE . IPLU) 
9  CONTINUE 

CALL  STEP(1PLB,1PULB,IPVLB, 

1  IPLA. IPULA, IPVLA, 

2  IPLA, IPULA, IPVLA, 

3  IP,IPU,IPV,IPUU,IPUV,IPVV, 

4  IPFX,IPFY,IPFXL,IPFTL,IPU,OT,F,CF,USX,USY) 

CALL  NAXNIN(IP) 

CALL  HAXMIN(IPU) 

CALL  NAXNIN(IPV) 

CALL  STEP(IPLC,IPULC,IPVLC, 

1  IPLA, IPULA, IPVLA, 

2  IPLB, IPULB, IPVLB. 

3  IP,IPU,IPV,IPUU,IPUV,IPVV, 

4  IPFX , IFFY , IPFXL . ZPFYL.IPW,DT, F,CF, USX, USY) 

CALL  COPY ( IPLA, IPLC ) 

CALL  COPY ( IPULA, IPULC ) 

CALL  COPY ( IPVLA , I PVLC ) 

10  CONTINUE 

IF  (PLOT)  CALL  0RIGIN(5.5,0. ) 

IF  (PLOT)  CALL  ENDPLT 

STOP 

END 

SUBROUTINE  INPUT ( STRING ,X  ) 

CHARACTER  STR ING <  6 ) 

UR  I TE ( OUTLUN , 1 )  STRING, X 

1  FORHAT <  '  Enter  value  for  '  ,6A1,  Perhaps' , 1 F 1 2. 3 ) 

ACCEPT  *,X 
RETURN 
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Connon  block  TRG.DCL 


COMMON  /TRG/  S' SMAX) , IV< IVMAX > . It < ICHAX > , 
.  SCALE, 

.  NEUV,NEUC,IVSIZ,ICSIZ,NUMC,NUHO,NBV,IBVF 


Library  of  subroutines  for  triangular  grid  codes 


f INCLUDE! 7HAC. MAC) (DEL TOK 

SUBROUTINE  ADVECT < IPCN, IPCL , IPC , IPU, IPU ,DT ) 
BINCLUDEITRG.DCL) ' 

J=1 

DO  1  I=1,NUNV 
K=IV( J) 

NV=IV( J+1 ) 

KIPC=K+IPC 

KIPU=K+IPU 

KIPU=K+IPU 

KIPCN=K+IPCN 

S(KIPCN)=S(KIPC)»(S(KIPU> *S(KIPU)fS<KIPUM >*S(KIPU+?>) 
DO  2  IN=1,NV 
JIN=J+ IN 
JP=IV(  JIN-H ) 

KP=I V I JP ) 

KIPCN=K+IPCN 
KPIPC=KP* IPC 
KPIPU=KP+IPU 
KIPUIN=K+IPU+IN 

2  S(KIPCN)=S(KIPCN)+S(KPIPC)*IS(KPIPU)»S(KIPMIN)* 

1  S<KPIPUH)»SIKIPUIN  +  7)) 

KIPCN-K+IPCN 

KIPCL=K+IPCL 

S!KIPCN)=S<K1PCL)-DT*S<KIPCN) 

I  J= J+IWSIZ 
RETUR’ 

END 


SUBROUTINE  CIRBND 
BINCLUDEITRG.DCL) ' 

DTH=6 . 2831 8A/NBV 

J*IBVF 

TH=0. 

DO  I  1=1, NBV 

IF  <J  .EO.  0)  GO  TO  9 

IF  < I V ( je 1 )  .EQ.  6)  GO  TO  ? 

K=IV< J) 

S(K)=COS<TH) 

SIKH  )=SIN(TH) 

TH=TH+DTH 
1  J=IV<  J  +  7 ) 

IF  (J  .EO.  0)  RETURN 

8  FORMAT  I "  CIRBND  ERROR',2110) 

9  UR1 TE ( OUTLUN , 8  )  I,J 
STOP 

END 
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SUBROUTINE  CUNF'L  I  ( IFLD, CONL'AL.N) 

(INCLUDE < TRG. DCL  )  “ 

DIMENSION  CONVAL (1),F(3),KU<3) ,XV<3) , YV(3) ,XL (2)  f U (2  ) 
J=1 

DO  1  1  =  1, NUNC 
DO  2  JJ=1 ,3 
JJJ=J+JJ 
JV=IC  <  J JJ ) 

KA=IV(JV) 

KV( JJ )=KA 
KAIFLD=KA+IFLD 
F ( JJ)=S!KAIFLD> 

XV< JJ)=S(KA) *SCALE 

2  TV < JJ)=S(KAH  )*SCALE 
DO  4  NC=1  ,N 

NI  =  0 

FC=CONVAL  <  NC ) 

DO  3  JA=1  ,3 
JB= JA* 1 

IF  <JB  .EQ.  4)  J8=1 

IF  (<F(JA)-FC)«<FC-F<JB>>  .LE.  0.)  GO  TO  3 
NI=NI+1 

ALPHA=(FC-F<  JA)  )/<F<  JB)-FUA1 ) 

XL (NI)=(1. -ALPHA MXV(JA)+ALPHA*XVUB) 

YL  <  NI )  =  <  1 . -ALPHA) *YV( JA ) ♦ALPHA*YV( JB ) 

3  CONTINUE 

IF  (NI  .EQ.  0)  GO  TO  4 
CALL  LINE<XL,YL,2,1,0,0) 

4  CONTINUE 

1  J= J+ICSIZ 
RETURN 
END 


SUBROUTINE  CONSRV ( IFROM) 

(INCLUDE ( TRG . DCL ) * 

SUN=0. 

J=1 

DO  1  1=1 ,NUNV 
K=IU< J) 

KIFROH=K*IFROH 

SUH»SUN+S(KIFROH) 

1  J= J+IVSIZ 

3  FORHAT ( "  FIELD' ,  15, '  INTEGRAL  , 1 F 15.P > 
URITE(0UTLUN,3)  IFRON,SUN 
RETURN 
END 
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SUBROUTINE  COPY ( I  TO, IFRQM) 
$  INCLUDE i  TRG . DCL ) ' 

J--1 

DO  I  1=1, NUHV 
K=IV( J) 

KITO=K+ITO 

KIFROH=K+IFROH 

3(KIT0)=S(KIFR0M) 

1  J=J+I VSI2 
RETURN 
END 


SUBROUTINE  CTU(IPOS) 

LOGICAL  NEIGH 
f INCLUDE ( TRG . DCL ) “ 

1  =  1 

DO  1  11=1 ,NUMV 
NV=IV( 1  +  1  ) 

K=IV( I ) 

X=S(K) 

Y  =  S  <  K  +  1 ) 

KIPOS=K+IPOS 

S(KIPOS)=0. 

DO  2  J=1,NV 
KIPOSJ=KIPOS+J 
2  S(KIP0SJ)=0. 

NVH1 =NV- 1 
DO  3  JP=1,NVH1 
I JP 1 = I ♦ JP+ 1 
MP= I V ( I JP 1 ) 

KP=IV(NP) 

PX=S(KP) 

PY=S(KP+1 ) 

JPP1=JP*1 
DO  3  JO=JPP1 , NV 
I JQ1 =1 ♦ JQ+ 1 
NQ  =  I V( I J  Q 1  ) 

IF  (.NOT.  NEIGH(NP,NQ> )  GO  TO  3 
KO=IV(NO) 

QX=S(KQ) 

QY=S(KQM ) 

PO»(PX-X)*<OY-T)-(PY-Y)=*(QX-X> 

S I  G=  1 . 

IF  (PO  .LT.  0.)  SIG*-1. 

IF  (IVUP+l)  . EQ.  6)  GO  TO  4 

DX».5*SIG*(PY-Y) 

DY=.3*SIG*(PX-X) 

KJP=KIPOS*JP 


■k _ _ _ _ _ _ _ JL 
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5  U'.IPOS )  =  S '  K1PUS  )  ♦ .  75* DX 
3tKlP0S  +  7)  =  S(KIP0S*7)  +  .75my 
S(KJP)-S(KJP)+.25*DX 
S<K JP*7 )=S(KJP+7) ♦,25*PY 

4  IF  (IV(JQH)  .£0.  6)  GO  TO  5 
DX--.5*SIG*(QY-Y) 

DY=.5*SIG*(QX-X) 
S(KIPQS)=S(K1P0S)+.75*DX 
S(KIP0S*7)=S(KIP0S+7) ♦ .75»DY 
K JP=KIPOS+ JO 
S<KJP)=S(KJPH.25*DX 
S(KJP*7)=S(KJP+7)+.25*DY 

5  CONTINUE 

DX--SlG*(FY-.5*(QY+Y))/3. 
tn=SIG*(PX~.5*<QX+X))/3. 
EX-SIG*(QY-.5*(PY*Y) )/3. 
EY=-SIG*(QX-.5*(PX»X))/3. 
S(KIP0S)=S(KIP0S)+j.*(DX+EX)/12. 
S(KIP0S  +  7)=S(KIP0S  +  7)t5.*(DY*-EY)/1 
KJP-KIPOS+ JP 

S(KJP)=S(KJP)*PX/6.+j.*EX/12. 

S(KJP+7)=S'KJP*7)+DY/6.+5.*EY/t2. 

KG=KIPOS+JO 

S ( KO ) =S ( KO) ♦  ). *DX/ 1 2. +EX/A. 
S(KQ+7)=S(KQ  •7)O.*0Y/12.+EY/A. 

3  CONTINUE 
I=I+1VSIZ 
1  CONTINUE 
RETURN 
END 


SUBROUTINE  CVLU(IPOS) 

LOGICAL  NEIGH 
♦INCLUDE) TRG.DCL) * 

A1=11./54. 

A2*?./108. 

1  =  1 

DO  1  11*1 ,NUHV 
NV=IV(I+1) 

K= I V < I ) 

X*S(K) 

r=s  <  k+ i ) 

KIPOS=K+IPOS 
S<KIPOS)=0- 
DO  2  J=1,NV 
KIPOSJ*K*IPOS+J 

2  S(KIP0SJ)=0. 

NVN1*NV-1 

DO  3  JP*1 rNVH1 
I JP=I+JP 
MP— I V  < I JP+1 ) 

KP=IV(NP) 

XP*S(KP) 

YP*S(KP+i  ) 

JPP1*JP+1 
DO  3  Jfl*JPP1,NV 
IJQ*I+JO 
NQ*IV( I J  Q  + 1 ) 

IF  (.NOT.  NEIGHINP, NO) )  GO  TO  3 
KQ*IV(NQ) 

XQ*S(KQ> 

YQ*S(KQ+1 ) 

DELTA*. 3*ABS( (XP-X) * ( YQ-Y >-< TP-Y) • (XQ-X) > 
KIPOS*K+IPQS 

S(KIPOS)=S(KIPOS)+A1  »DELTA 
KIPOSJP«KIPQS*JP 
S(K1P0SJP)*S(KIP0SJP)*A2»DELTA 
KIPOSJQ«KIPOS+JQ 
S  <KIPOSJO)*S<  KIPOSJQ  >  ♦A2*;PEL  I A 

3  CONTINUE 

I  *  I ♦ I VS  1 2 
1  CONTINUE 
RETURN 
END 
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SUBROUT INE  DRAUV< IPL, IPUL , IPVL f IP, If'U, IPV, ASCALE, 1PU ) 
*  INCLUDE <  TRG . DCL ) ' 

D1AENSION  X15>,Y<5) 

ERR=t.£-T 

CALL  UNLUNPT IPU , IPL , IP ,ERR ) 

CALL  UHLUHPUPU,  I  PUL  t I PU, ERR) 

CALL  UNLUMPC IPU, IPVL, IPV, ERR) 

J=t 

DO  10  1*1 ,NUHV 
K=IV( J) 

XX=S<K)»SCALE 
YY=S(KH  ) ♦SCALE 
KIP*K*1P 
P=S(KIP ) 

U-S(KIPU)/P*SCAL£ 

KIPV=K+IPV 

V*S(KIPV)/P*SCALE 

UV=SQRT<U*U+V*V> 

IE  <UV  .IT.  .01*ASCALE>  GO  TO  10 

DX=ASCALE*U*.5 

DY=ASCALE*V*.5 

X  ( 1  )  -  X  X  ♦  D  X 

X(2)=XX-DX 

Y(1)*YY+5Y 

Y<2)=YY-DY 

DX-. 4*0X 

DY=.4*DY 

X<3)=X<1)~DX-DY 

Y (3) *Y ( 1 ) -DY +DX 

X<5)*X( 1 >-DX*DY 

Y(5  )*  YU ) -DY-DX 

Y ( 4 ) *Y  ( 1 ) 

X(4)=X<U 

CALL  LINE<X,Y,2, 1,0,0) 

CALL  LIHE( X ( 3 ) , Y ( 3 ) , 3, 1,0,0) 

10  J=JMVSIZ 
RETURN 
END 


SUBROUTINE  DRCELL 
♦  INCLUDE ( TRG . DCL )  “ 

DIMENSION  X  <  2 ) ,  Y  ( 2 ) 

J=1 

DO  1  1=1 ,NUMV 
K  =  I V  (  J) 

XU  )  =  S  <  K  '  *SCAL£ 

TU  )=S( K+1  )*SCAL£ 

CALL  SYMBOL  < X , Y , .05,2,0. 

1  J= J+I VSI Z 
J=1 

DO  2  1=1, NUMC 
JA=IC( J  +  1  ) 

JB= IC ( J+3 ) 

JC=IC( J+3) 

KA=IV( JA) 

KB=IV(JB) 

KC=IV( JC) 

XA=S(KA) +SCALE 
YA=S(KAU  )*SCALE 
XB=S ( KB ) *SCALE 
YB=S ( KB+ 1 ) *SCALE 
XC=S(KC)*SCALE 
YC  =  S(KCU )*SCALE 
XG= ( XA+XB+XC ) /3  . 

Y  G  = ( YA+YB+YC ) /3. 

XU  )=XG 
YU  )  =YG 

X(2)=.5*(XA+XB) 

Y ( 2 )  =  .5* ( YA  +  YB ) 

CALL  LINE (X,Y, 2,1 ,0,0) 

X ( 2 ) = . 5* ( XB+XC ) 

Y ( 2 )  =  .5* ( YB  +  YC ) 

CALL  LINE(X,Y,2,1 ,0,0) 

X ( 2 ) = . 5*( XC+XA ) 

Y ( 2 )  =  .5*<  YX  +  YA ) 

CALL  LINE(X,Y,2,1 ,0,0) 

2  J=J+ICSIZ 
RETURN 
END 
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SUBROUTINE  DRGFflD 
$1NCLUDE(TRG.DCL) ' 

DIMENSION  X(2),Y(2) 

J=1 

DO  1  1=1 ,NUMV 
K=IV( J) 

NV=IV<  J+N 
X(1 )=S(K)*SCALE 
Y ( 1 )  =  S ( C+  1 )*SCALE 
DO  2  IP=1 ,NV 
JIP=J+ IP 
JP  =  IV<  JIP  +  1 ) 

IF  <JP  .GT.  J)  GO  TO  2 
KP=IV< JP) 

X ( 2 ) =S ( KP ) *  SC AL  E 
Y (2)=S(KP+1 )*SCALE 
CALL  LINE < X , Y,2, 1 ,0,0) 
2  CONTINUE 
1  J=J+IVSIZ 
RETURN 
END 


SUBROUTINE  EUPLOT  <  I^ZMAX  ,ZMIN) 

♦  INCLUDE (PAGE . DCL ) * 

LOGICAL  FAIL 

CHARACTER  ISP, ISTAR, IVERT, IHOR 

DATA  ISP/ '  '/, ISTAR/'*'/, IVERT/'! '/,IHOR/'_V 

DO  1  J=1 , MAPHAXV 

HO  2  1=1 ,MAPMAXH 

IMAP ( I , J ) = ISP 

IF  (I  .Ed.  1  .OR.  I  .EQ.  HAPNAXH)  WAP< I , J)=IVERT 
IF  (J  .EQ.  1  .OR.  J  .EQ.  MAPMAXV )  IMAP( I , J)=IHOR 

2  CONTINUE 
1  CONTINUE 

A=2./(MAPMAXH-1 .) 

B=- 1 . -A 

C= ( NAPHAXV-1 . )/(ZMIN-ZMAX) 

D= 1 . -C*ZNAX 
DO  3  1=1 ,MAPMAXH 
X=A*I+B 

CALL  INTP0L<IP,X,0.,Z,FAIL) 

IF  (FAIL)  GO  TO  3 
J=C*Z*D 

IF  (J  .LT.  1  .OR.  J  .GT.  HAPMAXV)  GO  TO  3 
IMAP ( I , J )  =  IST  AR 

3  CONTINUE 

4  FORMAT(')  ' ,/,  (5X,MAPHAXH  AD) 

URI TE ( OUTLUN, 4 )  IMAP 

RETURN 

END 


SUBROUT INE  FLLWV ( IPOSA ,  IPOSB , IPOSC ) 
(INCLUDE  (  TRG. BCL )  * 

J=1 

DO  1  1=1, NUHV 
K=IV( J) 

KIPQSA=K+ IPOSA 
KIP0S6=K+IP0?B 
K I P  0  S  C  =  K * IPOS' 

S(KIP0SA)  =  .5*(S( K IPOSB ) ♦S(KIPOSC) ) 

1  J=J>IUS1Z 
RETURN 
END 


SUBROUTINE  FRF <IPL , I  PUL , IPVL , IP , IPU, IPV, IPFX , IPF Y , IPF XL, IPFYL , 
.  DT,F,CF,USX,USY ) 

(INCLUDE ( TRG. DCL ) ' 

CALL  UNLUNP(2,IPL,IP,1 .E-2) 

CALL  UNLUrtP (2,IPUL,IPU,1.E_2) 

CALL  UNLUtIP  (  2  ,  IPUL ,  IPU ,  1 . E - 2 ) 

J=1 

DO  I  1=1 ,NUNM 

IF  (IV(J+1)  .NE.  6)  GO  TO  1 
K  =  1VU> 

KIP=K+IP 

P=S(KIP) 

KIPU=K+IPU 
U=S(KIPU)/P 
KIPV=K+IPU 
V=S(KIPV)/P 
UV=SQRT ( U*U+V*V) 

KIPFX=K+IPFX 

S(KIPFX)=(USX-CF»U>UV)"P 
KIPF  T=K+IPFY 

S(KIPFY)  =  (USY-CF*,/«UV)*P 
t  J-J+IVSIZ 

CALL  LUHP(2,IPFX,IPFXL) 

CALL  LUMP ( 2 , 1 PF  r . I PF YL  > 

J=1 

DO  2  1=1 ,NUHV 
K  =  I V  <  J  > 

A1 1  =  1 . 

A1 2=-DT  #F 
A21 =DT*F 
A22= 1 . 

KIPFXL=K+IPFXL 

KIPUL=K+IPUL 

B1=DT*S(KIPFXL)+S(KIPUL> 

KIPFYL=K+IPFYL 

KIPVL=K+IPVL 
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B2=PT*S(KIPFYL)+S<KIPVO 
D-fll 1  *  A  2  2 - A 1 1 *A21 
S<KIPUL)  =  <BUA22-A12*B2)/D 
S  <  K IPVL )  =  ( A1 1*B2-B1*A21  )/D 
2  J=J+IVSIZ 
RETURN 
END 


SUBROUTINE  GHC<IP) 

♦  INCLUDE ( TRG . DCL  > 

IC ( NEUC+1 ) = I P 

I C  <  NEUC>2 )  =  1F+ IVSIZ 

IC(NEUC+3)=NEUV-IVSIZ 

IP=IP+IVSIZ 

NEUC=NEUC+ICSIZ 

RETURN 

END 


SUBROUTINE  GUC(IP) 
♦INCLUDE! TRG. PEL  ) ' 

IC(NEUC+1 )=IP 

IC<NEUC+2)=NEWV 

IC!NEUC+3)=NEUV-IVS1Z 

NEHV=NEUV-*IVSIZ 

NEUC-NEUC+ICSIZ 

RETURN 

END 
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SUBROUTINE  HEXTOP(N) 

♦INCLUDE ( TRG.DCL) 

SCALE=4. 

DO  IOO  1=1 , ICHAX 

100  IC(I)=0 

DO  101  1=1 ,IVMAX 

101  IV(I)=0 
NEUC= 1 
!VSIZ=8 
I C  S I Z  =  4 
IP=1 

NEUV  =  1  +  01+2  )*IVSIZ 
IU=N 

DO  1  1=1, N 
DO  2  J=1,IU 
CALL  GUCUP) 

2  CALL  GDC < IP ) 

CALL  GUC(IP) 

IP=IP+IVSIZ 

NEUV=NEUV+IVSIZ 

1  IU=IU+1 
IU=IU-1 
DO  3  1  =  1, N 
DO  4  J=1,IU 
CALL  GDC ( IP) 

A  CALL  GUCUP) 

CALL  GDC ( IP ) 

IP=IP+IVSIZ 

NEUV=NEUV+IVSIZ 

3  IU=IU-1 
NEUV=NEUV-I VSIZ 
NUHC= (NEUC-1 ) /ICSIZ 
NUNVUNEUV-1 ) /IVSIZ 
WRITE ( OUTLUN , 5 )  NUHC,NUHV 

5  FORMAT < '  Allocated  cells  and  vertices' ,2110) 
JC=1 

DO  10  1=1, NUNC 
IVA=IC( JC+1 ) 

IVB=IC( JC+2) 

IVC=IC< JC+3) 

CALL  VLINK1 I V A , IV8) 

CALL  VL  INK< I VB, I VC ) 

CALL  VL1NK1 1VC,IVA) 

CAU  VLINK <  I VB,  IVA) 

CALL  VLINK<IVC,IVB) 

CALL  VLINK(IVA,IVC) 

10  JC= JC+ICSIZ 
IP*1 

DO  20  1*1 , NUNV 

IF  (IVUP+I)  .LT.  6)  60  TO  21 
20  IP*IP+IVSIZ 


J 


58 


22  FORMAT ( '  HO  BOUNDARY  POINTS' ) 

URI IE(OUILUN,22) 

STOP 
21  I B  = I P 
IBVF=IB 
IP=0 
NBV=1 

2?  IV1 IB+7)=1P 
NV=IV( IB  +  1 ) 

DO  23  1=1 rNV 
IBI 1 =IB+I+t 
J=IV(IBI1 ) 

IF  (J  .£0.  IP)  GO  10  23 
IF  U  .EQ.  IBVF )  GO  TO  25 
IF  (IVIJ+1)  .LT.  6)  GO  TO  24 

23  CONTINUE 

26  FORMAT ( '  BOUNDARY  ERROR  ) 

URITEtOUTLUN, 26) 

CALL  VCDUHP 
STOP 

24  NBV=NBV+1 
I  P=  IB 
IB=J 

GO  TO  27 

30  FORMAT < '  NUMBER  OF  BOUNDARY  VERTICES' ,110) 

25  URITE(0UTLUN,30)  NBV 
IBVF*IB 

RETURN 

END 


SUBROUTINE  INITF ( IP, IPU , IPV ) 
f INCLUDE! TRG.DCL  ) ' 

J*1 

DO  1  I = I , NUMV 
K  =  IV<  J) 

KIP*K+IP 
S(KIP)  =  t . 

KIPU*K*IPU 

S<KIPU>=0. 

KIPV*K+IPV 
S  <  K IPV  >  =  0  - 
1  J= J+IVSIZ 
RETURN 
END 


4 
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SUBROUTINE  INTPOU IFLD,X , Y , Z . FAIL > 

♦INCLUDE ( TRG . DCL ) ' 

LOGICAL  FAIL 
FAIL=. TRUE. 

I C  E  L  L  =  1 

DO  10  1=1 ,NUHC 
IP=IC( ICELL+1 ) 

IQ=IC ( ICELL+2 ) 

IR=IC ( ICELL+3 J 
JP  IV(IP) 

JQ=1V(IQ) 

JR= I V ( IR ) 

XP=S ( JP ) -X 

YP=S(JP+1 )-Y 

XQ=S ( JQ ) -X 

YQ=S ( JQ+ 1 )-Y 

XR=S(JR)-X 

YR=S( JR+1 )-Y 

PQ=XP*YO-YP*XQ 

QR=XQ*YR-YQ*XR 

RP=XR*YP-YR*XF 

IF  (PQ* OR  .LT.  0.)  GB  TO  t 

IF  (QR*RP  .LT.  0.)  GO  TO  1 

IF  <PQ*RP  .LT.  0.)  GO  TO  1 

JPIFLB=JP+IFLD 

ZP=S( JPIFLD) 

JQIFLD= JQ+IFLD 
ZO=S< JQIFLD) 

JRIFLD=JR+IFLD 
ZR=S( JRIFLD ) 

D=XP*(YQ-YR )  -YPM  XQ-XR ) ♦XQ*YR-YG*XR 

Z  =  <XP#<YQ*ZR-ZO*YR)-YP*<XO*ZR-ZO=»XR)+ZP»(XQ*YR-YO*XR))/Il 
FAIL=. FALSE. 

RETUFL 

1  ICELL=ICELl+ICSIZ 
10  CONTINUE 
RETURN 
END 


t 


i 
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FUNCTION  INTRI <  X  f  Y  , 1  CELL  ) 

♦  INCLUDE ( TRG . DC L ) A 
LOGICAL  INTRI 
I P= I C ( I  CEL  L ♦ 1 ) 

IQ=IC  ( ICEI.L  +  2 ) 

IR=IC( ICELL+3) 

JP=IV( IP ) 

JO=IV( IQ) 

JR=1V<  JR) 

XP=S(JP)-X 

YP=S( JP  +  1 ) -Y 

XO=S(JO)-X 

YR*S( JQ+ 1 ) -Y 

XR=S( JR)-X 

YR=S( JR+1  )-Y 

PQ=XP*YQ-YP*XO 

QR=XQ*YR-YQ*XR 

RP=XR*YP-YR*XP 

IF  (PO*QR  .LT.  0.)  GO  TO  1 

IF  (QR*RP  .LT.  0.)  GO  TO  1 

IF  <PQ*RP  .LT.  0.)  GO  TO  1 

I HTR I = . TRUE. 

RETURN 

1  INTRI=. FALSE. 

RETURN 

END 
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SUBROUTINE  LUMP < IPOS, IFRON, ITO) 
♦INCLUDE! TR6.DCL ) * 

J  =  1 

DO  1  1=1 ,NUMV 
K=IV< J) 

NV=IV<  J  +  1 ) 

KITO=K*ITO 
KIP0S=K*IPOS 
KIFROM=K+IFRON 
S  ( K I  TO )  =S  <  KIPOS  )^S  <  KXFROli) 

DO  2  IN= I ,NV 
JINl=J+IN+t 
JP=IV( JIN1 ) 

KP=IV( JP) 

KIPIN=KIPOS+IN 

KPIFR=KP+IFROM 

2  S<KITO)=S(KI TO )+S (KIRIN )»S(KPIFR) 
1  J=J*IVSIZ 
RETURN 
END 


SUBROUTINE  NAXNIN(IFLD) 

♦  INCLUDE! TRG.DCL)  ‘ 

XNAX=S< IFLD+1 ) 

XNIN=XNAX 

J=1 

DO  1  1=1 ,NUHV 
K=IV( J) 

KIFLD=K+IFLB 
X=S  <  K IFLD  > 

IF  (X  .GT.  XNAX )  XNAX=X 
IF  (X  .LT.  XNIN)  XNIN=X 

1  J= J+IVSIZ 

2  FORMAT (  '  MAX, MIN''  , 2F1 0 .6 ) 
URITE(0UTLUN,2)XMAXfXHIN 
RETURN 

END 


FUNCTION  NEIGH( IP, IQ) 
*INCLUDE<  TRG.DCL ) * 

LOGICAL  NEIGH 
NV=IV( I P+ 1 ) 

DO  1  1=1,  NV 
IPI1=IP*I+1 

IF  (IV(IPII)  ,EQ.  IQ)  GO  TO  2 

1  CONTINUE 
NEIGH=. FALSE. 

RETURN 

2  NEIGH=.TRUE. 

RETURN 

END 


SUBROUTINE  PRCON(IFLD,XHIN,XMAX> 
LOGICAL  FAIL 

DINENSION  ICON < 13> ,LIN( NAPHAXH) 
DATA  1COH/'-' ,  0  ,'1' ,'2'  ,'3' ,'A' 
BX=-(HAPHAXH+1 . )/HAPNAXH 
BY=(HAPHAXV+1 . )/HAPHAXV 
URI TE ( OUTLUN , 5 ) 

DO  1  J= 1 , HAPHAXV 
Y=-J/(HAPHAXV/2. ) *B  Y 
DO  2  1=1 , NAPHAXH 
X=I/(HAPHAXH/2.)+BX 
R=X*X+Y*Y 
K=1 3 

IF  (R  .GT.  1.)  GO  TO  3 
CALL  INTPOL(IFLD,X,Y,Z,FAIL) 

IF  (FAIL)  GO  TO  3 
K=10.*(Z-XNIN)/(XHAX-XHIN) 

IF  (K  .LT.  0)  K  =  -1 
IF  (K  .GT.  10)  K=1 0 
K=K+2 

3  LIN< I )  =  I  CON ( K  ) 

2  CONTINUE 

4  FORHATOX, NAPNAXH  A1) 

URI TE ( OUTLUN, 4 )  LIN 

1  CONTINUE 

5  FORMAT ( "  1  ' ) 

RETURN 

END 
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SUBROUTINE  PRFLDS( It ,12) 

♦  INCLUDE ( TRG . DCL )  * 

DIMENSION  VUO) 

UR  I T  E  < OUTLUN  T 4 )  11,12 
4  FORHAT < '  F IELDS " f 1 4 , "  THROUGH'  , 14  > 
URITE(0UTLUN,4)  11,12 
J=1 

DO  1  1  =  1,  NUHV 
K = I V  ( J ) 

NN=1 

DO  2  N=1 1 , 12 

KN=K+N 

V(NN)=S(KN) 

NN=NN+1 

IF  (NN  .LE.  10)  GO  TO  2 
UR  I TE  <  OUTLUN , 3 )  I,V 
3  FORMAT ( '  ' , 18, 1  OF  1 2.6) 

NN=1 

2  CONTINUE 
NN=NN-1 

IF  (NN  .GT.  0)  URITE(0UTLUN,3)  I, (V(N),N=1 ,NN) 
1  J=J+IVSIZ 
RETURN 
END 


SUBROUTINE  QPRCON(IFLD,XMIN,XMAX) 

CHARACTER  IC0N(13) 

♦  INCLUDE ( TRG . DCL  )  ~ 

I INCLUDE (PAGE. DCL)' 

DATA  ICON/'-VOVI  V2V3VAV5V6V7V8V?' 

DO  2  J=1 , HAPHAXV 
DO  2  1=1, NAPNAXH 
2  I  NAP ( I , J )  =  I  CON ( 1 3 ) 

IVERT=1 

DO  1  I 1=1 ,  NUHV 
KK=IV( IVERT ) 

X=S(KK ) 

Y  =  S(KK  +  1  ) 

KKIFLD=KK+IFLD 

Z=S(KKIFLD) 

1=  (HAPHAXH/2.)«X*HAPHAXH/2.M. 

IF  (I  .LT.  1)  1=1 

IF  (I  .GT.  HAPHAXH)  I=HAPHAXH 

J*-(HAPHAXV/2. >*Y*HAPHAXV/2.+1 . 

IF  (J  .LT.  1)  J=1 

IF  (J  .GT.  HAPHAXV)  J=HAPHAXV 

K* 1 0. *<  Z-XHIN)/( XHAX-XHIN) 

IF  (K  .LT.  0)  K=- 1 
IF  <K  .GT.  10)  K»10 


'/ 


i 

>i 
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K=K+2 

IMAP  < I, J)=ICON(K) 

I  I VERT=I VERT+IVSIZ 
3  FORMAT  ( ' 1 ' ,/,  (5X,HAPMAXH  AD) 
URITE(OUTLUN,3)  IMAP 
RETURN 
END 


SUBROUTINE  RIGBNDdPUL,  IPVL, IPU , IPV) 
$  INCLUDE ( TRG . DCL  ) " 

CALL  UNLUNP(2,IPUL,IPU,1.E-5> 

CALL  UNLUMPI2, IPVL, IPV, 1.E-5) 

J=IBVF 

DO  1  I = 1 , NBV 
K=IV( J) 

KIPU=K*IPU 

S<KIPU)=0. 

KIPV=K+IPV 

S(KIPV)=0. 

1  J=IV( J  +  7) 

CALL  LUMP! 2, IPU, IPUL ) 

CALL  LUHP<2, IPV, IPVL) 

RETURN 

END 


SUBROUTINE  SBROT ( ONEGA, IPOSU ) 
I INCLUDE! TRG. DCL )  * 

1  =  1 

DO  1  11=1, NUMV 
K=IV( I ) 

X=S(K) 

r=s<K+i) 

I=I*IVSIZ 
KIP0SU=K+IP0SU 
S(KIPOSU)*~Y*OHEGA 
1  S ( K I POSU ♦ 1 )*X*ONEGA 
RETURN 
END 
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SUBROUTINE  STEPTIPLN, IPULN, IPVLN , IPLO, IPULO. IPVLO. IPL. IPUL , IPVL, 
.  IP,IPU,IPV,IPUU,IPUV,IPVV,IPFX,IPFY,IPFXL,iPFYL,lPU,DT, 

.  F,CF,USX,USY> 

♦INCLUDE (TRG.DCL) 

CALL  UNLUNP(2, IPL,IP,1.E~2) 

CALL  UNLU«P<2,IPUL,IPU,1.E-2) 

CALL  UNLUNP(2, IPVL, IPV, 1 .E-2) 

J=1 

DO  1  1=1 ,NUNV 
K=IV( J) 

KIP=K+IP 

KIPU=K+IPU 

KIPV=K+IPV 

KIPUU=K+IPUU 

KIPUV=K+IPUV 

KIPVV=K+IPVV 

P=S(KIP ) 

U=S<KIPU)/P 

V=S(KIPV)/P 

S(KIPUU)=P*(U*U+.5*P) 

S<KIPUV)=P*U*V 

S<KIPVV)=P*<V*V*.5*P> 

1  J=J+IVSIZ 

CALL  TRNSPT < IPLN, IPLO , IPU , IPV , IPU, DT ) 

CALL  TRNSPT (IPULN, IPULO, IPUU.IPUV, IPU, DT) 

CALL  TRNSPT < IPVLN, IPVLO, IPUV, IPVV, IPU, DT ) 

CALL  FRF ( IPLN , IPULN, IPVLN, IP , IPU, IPV, IPFX, IPFY , IPFXL, IPFYL,DT , 

.  F,CF,USX,USY) 

CALL  R I GBND< IPULN, IPVLN, IPU, IPV) 

RETURN 

END 


SUBROUTINE  SUGRID1 NITER) 

♦  INCLUDE*  TRG. BCL ) ' 

RELPA= 1 . 388 
DO  9  MI  =  1  ,NI TER 
J=1 

CHMAX=0 . 

DO  1  1=1 ,NUHV 
NV=IV( J+l ) 

IF  ( NV  .LT.  6)  GO  TO  t 
X=0. 

Y=0. 

DO  2  K=1 ,6 
JK1 = J+l +K 
JN=IV( JK1 ) 

JS=1V(JN) 

X=X+S(JS) 

2  r=r+s<js*n 

JS=IV( J) 

XO=S< JS) 

ro=s(js+t ) 

X=X/6 . 

r=r/6. 

CH=SORT(  <X-X0):**2*(Y-Y0):**2) 

IF  < CH  .GT.  CHHAX )  CHHAX=CH 
SUS)=RELPA*X*<1  ,-RELPA)*SUS) 

S< JS-M  )=RELPA*YM1-RELPA)*S<JS+1  > 

1  J=J+IUSIZ 

8  FORMAT < '  MAXIMUM  DISPLACEMENT' , 1 F12. 8) 
URI TE ( OUTLUN , 8 )  CHHAX 

9  CONTINUE 
RETURN 
END 
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5UBF.0U1  INE  TRNSF'I  ( INF  U .  IOLD,  II  X,  IF  t ,  IPU,  DT ) 

♦INCLUDE ( TRG . DCL ) 

J=l 

DO  t  1=1 ,NUHV 
K  =  I V  < J) 

NV=1V(J*1) 

KINEU=K+INEU 
KIFX-K+IFX 
KIPU=K+IPU 
KIF  Y  =  K*  IFY 

S(KINEU)=S(KIFX)*S(KIPU)+S(KIFY)*S(KIPU+7) 

DO  2  IN=1 ,NV 
JIN= J+ I N 
JF'=IV(  JIN+1  ) 

KP= I V  <  JP  » 

KPIFX=KP+ IFX 
KIPUIN=KIF'U  +  1N 
FYPIF  Y=KP  +  IFY 

2  S<KINEU)=S(KINEU)+S<KPIFX)*S<KIPUIN>+S(KPIFY) *S(KIPUIN+7> 
KIOLD=K+IOLD 

S(KINEU)=S(KIOLD)-DT*S(KINEU) 

I  J=J+IVSIZ 
RETURN 
END 


SUBROUTINE  TSTCTU C IPOS ) 

*INCLUDE<TRO.DCL) " 

URITE(OUTLUN,A) 

4  FORMAT (  '  TEST  T-WEIGHTSV8X , ' I '  ,8X,  NV'  r ) 8X , 'TEST  X',HX,'TEST  Y  ) 
1  =  1 

HO  1  11=1, NUMV 
N V= I V  < 1  +  1  ) 

K  =  I V  ( I ) 

KIPOS=K+IPOS 

X=S(KIPOS) 

Y=S(KIP0S+7> 

DO  2  J=1 ,NV 
KIPOSJ=KIPOS+ J 
X=X+S  <  K IPOS J  > 

Y=Y+S<KIP0SJ+7) 

2  CONTINUE 

3  FORMAT < 1 X, 21 10, 2F20. 1 0> 

URITE(0UTLUN,3)  I,NV,X,Y 

1  I  =  IUl>SIZ 
RETURN 
END 


be 


SUBKOUi  lNfc  UNLUrtl'-'l  II-OS .  I  FROrt ,  I  TO ,  ERR ) 
t INCLUDE  (  IKCi.DCL)  * 

RElPA^I  .3 
NITER-25 

DO  1 1  N I T  =  1 .NITER 
CHHAX=0. 

J=1 

DO  10  1=1, NUrtV 
K  =  I V  ( J ) 

KITO=K+ITO 
SOLD  =  S ( K I  TO ) 

K IFROM=K+ IFROM 
S ( K I  TO ) =S ( K IFROM ) 

NV= I V ( J  +  1  ) 

DO  12  I N= 1 ,NV 
JIN-J+IN 
JP  =  IV<  JIN+I  ) 

KP- 1 V ( JP ) 

KIPSIN=K+IPOS+IN 
KPI TO=KP+ITO 

12  S(KITO)=S(KITO)-S(KIPSIN)*S(KF'ITO' 
KIPOS=K+IPOS 

S(KITO)=S(KITO)/S(KIPOS) 
CH=ABS(S(KITO)-SOLD) 
S(KITQ)=RELPA*S(KIT0)+(1 . -RELPA)*SOLD 
IF  ( CHhAX  .LT.  CH)  CHHAX=CH 

10  J=J+IVSIZ 

IF  (CHHAX  .LT.  ERR)  GO  TO  H 

11  CONTINUE 

13  FORHAT (  MAXIHUH  CHANGE  IN  UNIUMP  ,1F12.8) 
WRITE ( OUTLUN, 1 3 >  CHHAX 

M  CONTINUE 
RETURN 
END 


SUBKOUI INL  VCDUMP 
I  INCLUDE ! TRG. DCL  ) 

UR  I TE ( OUTLUN , 1 ) 
t  FORMAT!  •'  CELL  LIST'' ) 

J=1 

DO  2  1=1 ,NUHC 
JL=J*ICSIZ-1 

WRITE! OUTLUN, 3)  I , J , ( IC ( K ) , K= J, JL ) 

2  J= J+ICSI Z 

3  FORMAT ! 1 X ,61 1 0 ) 

J=1 

A  FORMAT ( ■'VERTEX  LIST' ) 

WRITE! OUTLUN , A ) 

DO  5  1=1, NUMV 
JL= J+IVSIZ 

WRITE (OUTLUN, 6)  I , J, ( IV(K  > , K=J, JL ) 
6  FORMAT ( < 1 X , 10110)) 

5  J=J+IVSIZ 
RETURN 
END 


SUBROUTINE  VLINK! IVA, IVB) 
(INCLUDE ( TRG . DCL )  * 

NV=IV( IVA  +  1 ) 

IF  < NV  .EQ.  0)  GO  TO  1 
DO  2  1=1, NV 
IVAI1=IVA*I+1 

IF  (IV(IVAII)  .EO.  IVB)  RETURN 
2  CONTINUE 
1  IV! IVA  +  1 )=NV+1 
IVANV=I VA+NV 
IV! IVANV+2)=IVB 
RETURN 
END 


SUBROUTINE  VSTO! IORIG,NUDS) 
IINCLUDE(TRG.DCL)* 

K  = I  OR  I G 
J=l 

DO  1  1=1, NUMV 
IV! J)=K 
K=K+NWDS 
1  J=J+IVSIZ 
RETURN 
END 


END 

DATE 
FILMED 


DTIC 


